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We have developed a general methodology for numeri- 
cal computations of electromagnetic (EM) fields and forces 
in matter, based on solving the macroscopic Maxwell's equa- 
tions in real space and adopting the Maxwell Stress Tensor 
formalism. Our approach can be applied to both dielectric 
and metallic systems of frequency-dependent dielectric func- 
tion; as well as to objects of any size and geometrical proper- 
ties in principle. We are particularly interested in calculating 
forces on nanostructures. 

Our findings confirm that a body reacts to the EM field by 
minimising its energy, i. e. it is attracted (repelled) by regions 
of lower (higher) EM energy. When travelling waves (of real 
wavevector) are involved, forces can be additionally under- 
stood in terms of momentum exchange between the body and 
its environment. However when evanescent waves (of com- 
plex wavevector) dominate, the forces are complicated, often 
become attractive and cannot be explained by means of real 
momentum being exchanged. 

We have studied the EM forces induced by a laser beam on 
a crystal of dielectric spheres of GaP. We observe effects due 
to the lattice structure, as well as due to the single scattering 
from each sphere. In the former case the two main features are 
a maximum momentum exchange (and largest forces) when 
the frequency lies within a band gap; and a multitude of force 
orientations when the Bragg conditions for multiple outgoing 
waves are met. In the latter case the radiation couples to the 
EM eigenmodes of isolated spheres (Mie resonances) and very 
sharp attractive and repulsive forces occur. Depending on the 
intensity of the incident radiation these forces can overcome 
all other interactions present (gravitational, thermal and Van 
der Waals) and may provide the main mechanism for forma- 
tion of stable structures in colloidal systems. 

PACS numbers: 42.70Qs, 41.20.Bt, 42.25Fx, 82.70.Dd 



I. INTRODUCTION 



in homogeneous media. Another startling effect is the 
possibility for no allowed optical modes to exist at cer- 
tain frequency ranges (band gaps) which, coupled with 
the ability to inject optical modes at specific frequencies 
by introducing defects in the crystalline structure, has fu- 
elled prospects for novel optical phenomena such as pho- 
ton localisation, inhibition of spontaneous emission and 
enhanced nonlinear effects. Novel devices (single-mode 
LEDs, thresholdless lasers, etc.) are also expected. 

Where do we find photonic crystals ? Natural gem- 
stone opals owe their beautiful iridescent colours to their 
structure which strongly diffracts visible light. However 
their relatively large amount of disorder as well as the 
poor refractive index contrast they display renders them 
impractical for the detailed tailoring of optical properties 
one needs to attain if the above mentioned novel effects 
are to be sought after. The first man-made photonic crys- 
tals were made ten years ago by Yablonovitch and they 
operated at the scale of microwaves. Intensified efforts for 
the past years have born fruit and we are now starting 
to have artificial opals at visible wavelengths, fabricated 
by a wealth of experimental techniques which bears wit- 
ness to the ingenuity with which the problem has been 
tackled. Some methods aim at building crystals from col- 
loidal suspensions utilising the influence of gravitational, 
electromagnetic or biochemical interactions, while oth- 
ers concentrate on selectively removing material from a 
matrix (etching) and thus forming a desired structure. 
'Mixed' approaches are also used (removing material fol- 
lowed by ordering or stacking of crystal layers) ||]. 

The theoretical understanding of photonic structures 
has centred at solving Maxwell's equations for the EM 
wavefield in the presence of matter. Thus the classical 
theory of vector fields is applied. This methodology, ap- 
propriate for scales large enough compared to the under- 
lying atomic structure, has been successfully applied to 
predict both the dispersion relation (band structure) as 
well as the scattering properties of such systems. 



A. Photonic structures 



Photonic materials possess structure (crystalline or- 
der) at the scale of the wavelength of light which they 
therefore diffract strongly, just like x-rays are diffracted 
from atomic crystals. By adjusting the geometric and 
dielectric parametres of a photonic crystal, we can alter 
the optical properties at will in the same way as a 

semiconductor crystal can be built to support desirable 
electronic properties. For example, the dispersion rela- 
tion (band structure) of light inside a photonic crystal is 
very different from the linear u) = kc behaviour we have 



B. Overview of photonic forces 

Nature favours phenomena that are energetically prof- 
itable. For example, the Van der Waals force between 
two macroscopic objects can be found by calculating the 
change in the self-energy of the vacuum-fluctuation wave- 
fields which multiply scatter off the objects. In electronic 
band structure, total energy calculations combined with 
variational arguments have been successfully used to pre- 
dict the equilibrium structure of crystals and to realisti- 
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cally estimate a plethora of physical quantities |j] . But 
while the former example involves virtual light waves and 
the latter employs real electronic waves, there is an in- 
creasingly widening range of systems where real photons 
apply; of those, photonic crystals form a subset whose 
potential for new physics and new applications is only 
beginning to be realised. 

In such a crystal, the photonic band structure con- 
tains the energy levels accessible to photons. But unlike 
atomic crystals, where the sea of electrons occupies a 
range of modes, it is possible with an external source of 
monochromatic light to selectively populate a particular 
mode in a photonic crystal. The light source injects many 
photons but owing to the bosonic and non-interacting 
nature of light, each photon would remain in the same 
state. Thus photonic crystals are suitable for observing 
single-mode band-structure effects, a property implying 
rich force spectra. To give an example, consider what 
happens as we traverse with a laser of tunable frequency 
the first band gap in a photonic crystal made of dielectric 
spheres in air. On the low energy side, the EM wavefield 
is concentrated outside the spheres {i.e. at regions of 
low dielectric constant), hence the energy of this mode is 
lower. Inside the gap, the fields are everywhere zero since 
all light is expelled from the crystal. And on the high 
energy side, the field is concentrated inside the spheres. 
This shifting of the wavefield from one region to another 
translates into changes in the spatial distribution of EM 
energy inside the crystal, and would therefore strongly 
affect the forces. 

Apart from performing total energy calculations, EM 
forces can be also computed in an alternative, more for- 
mal manner: namely within the Maxwell Stress Tensor 
(MST) formalism. In this methodology, the various com- 
ponents of the MST are constructed from the electric and 
magnetic fields in real space. Integration of the tensor 
around a closed surface surrounding the body of inter- 
est yields the total force and torque acting on it. There 
is in principle no restriction as to the size and shape of 
the body, nor as to its dielectric properties. The general- 
ity of this method makes its applicability wide, ranging 
from purely electric to purely magnetic to fully electro- 
magnetic effects; from virtual to real photons; and from 
travelling to evanescent wavefields. 

Our approach in this paper is as follows: First we ex- 
plain in some detail how the concept of a stress tensor 
arises in the calculation of forces. We define the tensor for 
the case of electromagnetism and show how it can be in- 
corporated in numerical calculations of photonic effects. 
We then put our methodology at work. As a demon- 
stration of some principal features of EM forces we treat 
first the case of a beam of laser light scattering from 
an homogeneous medium (metallic or dielectric) of vari- 
able thickness. This is followed by the main body of our 
results which regards light-induced forces on dielectric 
crystals. 



II. METHODOLOGY 
A. Theoretical framework 

1. The need for a stress tensor 

We are interested in calculating EM forces on macro- 
scopic bodies in the nanometre scales. It is therefore 
natural to seek expressions for the forces in terms of 
macroscopic quantities (such as dielectric functions). We 
therefore make at the outset an implicit adoption of the 
continuum picture of matter. We follow the treatment of 
Landau and Lifshitz 

A force acting on a macroscopic body is a sum of the 
forces fi,(r) on each of the volume elements that consti- 
tute the body. Thus the total force F is 

F = J {,{r)dv. (1) 

If {y can be expressed as the divergence of a tensor of 
rank two, then we can transform the volume integral into 
a surface integral and the total force becomes 

■^"^/^^'^"' <^^f3 = {x,y,z} , (2) 
= jT^f3dSp, (3) 

where a sum over repeated indices is implied from here 
onwards. The stress tensor Ta/s contains all the informa- 
tion we need. Note it has dimensions of energy density. 
This is not accidental, since forces can be also thought 
to arise due to spatial variations in the energy (energy- 
gradient picture) . They act in order to lower the system's 
energy i. e. to push the system into regions where the en- 
ergy is lowest. Therefore only the spatially varying part 
of TaB contributes to forces, a fact which can be seen in 
Eq. (g). Eq. (||) is of great practical benefit in that it is 
often easier to calculate the stress tensor on the surface 
of a solid body rather than in its interior. 

2. The stress tensor for electromagnetism 

So far the discussion is general. To concentrate on EM 
forces we need the stress tensor for the case of electro- 
magnetism. When a solid body is immersed in a linear, 
isotropic medium of permittivity eoCm and permeability 
MoMm, then the spatially non-constant part of the stress 
tensor in the region external to the body is given in terms 
of the electric (E) and magnetic (H) fields by 

Tap = ^Q^mEa Ep ^ ^iQ^rnHa Hp 

— T^^afji^O^mEjE^+^OflniHjHj) (4) 
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provided the medium remains in mechanical and ther- 
mal equilibrium under the influence of the EM fields {i.e. 
electrostriction and magnetostriction effects are ignored). 

Thus by integrating Ta0 of Eq. (^) over a closed sur- 
face S enclosing the body, the full EM force on the body 
is obtained. No further approximations or limiting as- 
sumptions need be made, as has often been done with 
other methodologies employed for EM effects (e.g. dipole 
approximation, vanishing particle size compared to the 
wavelength of light in scattering phenomena, etc.). Our 
methodology can be applied to dielectric or metallic ob- 
jects of any shape and at any frequency range in principle. 
The only requirement for meaningful results is that the 
integration is carried out over scales which are large with 
respect to atomic sizes, so that the continuum picture of 
matter is justified. For nanostructure scales of interest to 
us (> 50 nm) this requirement is surely satisfied: the dy- 
namics of the EM wavefield is more properly described by 
the equations of Maxwell rather than by the Schrodinger 
equation for all the electronic charges which make up the 
dielectric response of the system. Finally, our choice of 
macroscopic scales dictates a relevant timescale also (no 
macroscopic object can respond synchronously to THz 
oscillations): we must time- average the forces {i.e. the 
stress tensor) over many cycles when dealing with har- 
monic fields at optical frequencies. 



H(r,t) = HL E (^i^M + ^k'.') ^k'.' e^I'-'-'— '1 . 

fe',C7' 

(5b) 

In Eq. (^ E"„^,H"„^ are the amplitudes of the incident 
fields and e^'^./ are the unit vectors for each k' and each 
polarisation a' . The Kronecker delta functions serve as 
to add the contribution from the incident field, which is 
present on the —z side, to the reflected field. For the 
total fields on the +z side the only contribution comes 
from the transmitted fields. 

Having constructed the total fields outside the crystal, 
we invoke once again the Matrix Transfer equations to 
find the total field at each mesh point inside the struc- 
ture. When this is accomplished, the stress tensor com- 
ponents are calculated and time-averaged. The latter 
procedure can be made in a standard way since the fields 
are expressed as sums of plane waves at the same fre- 
quency uj. Finally, the forces are found by integrating 
the stress tensor over a surface enclosing the body of 
interest, which can be of arbitrary shape and complex- 
ity. In this work we present force calculations on isolated 
spheres, crystal layers and the entire crystal sample. 

III. RESULTS AND DISCUSSION 



B. Numerical methods 

We discretise Maxwell's equations on a simple cubic 
(SC) mesh. For calculations of photonic band structure, 
periodic boundary conditions are used in all directions. 
For transmission results, we consider a crystal finite in 
the z direction, having a thickness of N unit cells (lay- 
ers), N being allowed to vary. We shall discuss here the 
case where a monochromatic light beam of frequency oj, 
polarisation a and wavevector k is incident on the crys- 
tal from the —z side. Using the Transfer Matrix formal- 
ism we calculate reflection coefficients foi' scattering 
into plane waves of wavevector k' = ky -I- gy -I- zfc^ and 
polarisation cr', where k|| is the component of k parallel 
to the xy plane and g|| is a two-dimensional reciprocal 
lattice vector also lying on the xy plane. Our computer 
codes have been published and issues of stability and ac- 
curacy have been discussed elsewhere |l0| . 

The new ingredient is the calculation of forces. From 
the reflection coefficients we construct the total fields on 
the —z side of the crystal, 

E(r,t) = + Rk'.') e.v e^^'"—*! , 

k',<j' 

(5a) 



A. Homogeneous system - travelling waves 

Probably the simplest system for studying macroscopic 
radiation forces on materials is a homogeneous half-space 
(dielectric or metallic) upon which a single beam of light 
is incident. The simplicity of this model system will allow 
us to put forward in a transparent manner some impor- 
tant features of our methodology. Furthermore, as this 
problem is analytically solvable it provides a nice test- 
ground for our numerical results. 

We model the semi-infinite medium by considering a 
layer of material of sufficient thickness in the z direc- 
tion so that all light entering the layer is absorbed before 
reaching the far end. For a metal this would only take 
a thickness equal to the skin depth (a few nm typically) 
when Lo lies below the plasma frequency ujp; for a real 
dielectric, the tiny absorption present makes it necessary 
to consider much thicker samples, perhaps as much as 
several cm. In either case, only the incident and specu- 
larly refiected wavefields are present outside the system. 
The force on the layer is found by integrating the stress 
tensor for the total field (incident plus reflected) on a sur- 
face Sz which lies on the —z side immediately above the 
medium and extends to infinity in the xy plane (the con- 
tribution from the integration on the far end is zero, since 
the fields vanish there). A more appropriate quantity for 

this system is the EM pressure Pa — {1/Sz) § TazdSz on 
the layer. 

As a dielectric system we chose GaP which is transpar- 
ent {i.e. has low absorption and the dielectric function e 
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is mostly real and positive) over near infrared and visible 
frequencies but behaves as a poor metal {i.e. there is a 
plasma frequency at ~ 13.5 eV but the absorptive part 
of e becomes significant with respect to the real part) in 
the ultraviolet spectrum. As a metallic system we use 
Al, for which the dielectric function can be modelled as 

UJ?. 

e = 1 - , V , , (6) 
uj{lu + ?7) 

where hujp — 15 eV and the damping coefhcicnt for the 
intraband transitions is = 0.1 cV. Our data for the 
dielectric properties of these materials come from |l^ . 

We summarise our main findings for the forces on this 
system: 

(a) A homogeneous system can only scatter light into 
directions belonging to the scattering plane (the plane 
defined by k and z) so EM forces must lie on this plane. 
Also the polarisation of light is conserved upon scatter- 
ing. 

(b) The response to light incident on a homogeneous half- 
space will range from full reflection (e.g. for a perfect 
metal at frequencies below ujp) to full absorption (for a 
dielectric at the Brewster angle and p-polarised light). 
For normal incidence, the pressure on the half-space will 
range respectively from 2Pi„c (full reflection) to Pine (full 
absorption), where Pine — Io/cq is the radiation pressure 
of the incident beam which has intensity Iq and travels 
at the free-space speed cq. If we replace the half-space 
with a layer of finite thickness such that some light is 
transmitted through the far end, taking care to include 
the contribution from the far end in the stress tensor inte- 
gration, then the lower bound for the force becomes zero, 
corresponding to full transmission (it i.e. zero reflection 
and zero absorption) . These effects are demonstrated for 
GaP and Al samples in Figs. |l|, ||. It is clear from these 
plots that the force on the system equals minus the net 
momentum change of the light, just as one would expect 
from the action-reaction principle. Thus the explanation 
of forces in terms of the momentum-balance picture is 
valid. 

The three values {0,Pinc and 2Pinc) constitute the nat- 
ural limits for total {i.e. external) EM forces on a macro- 
scopic body (homogeneous or not, finite or not) of ar- 
bitrary size, due to incidence of a single beam of light. 
The only way for these limits to be exceeded is when a 
resonance condition is fulfilled. We shall come to this 
interesting case when dealing with forces on dielectric 
crystals, (c) Since forces can just as well be attributed 
to spatial variations of energy as to temporal variations 
of momentum, we can use the energy-gradient picture to 
study EM forces on bodies. In this picture, the EM force 
on an infinitesimal volume element located at r is 



of material extending from to d in the z direction, the 
total {i.e. average) force on the layer is 



-VU . 



(7) 



where U is the energy of the EM wavefield at r. The 
force on a macroscopic body is found by adding the con- 
tributions from all elements comprising the body. For ex- 
ample, when light is incident at normal angle on a layer 



V,U dz = - z ^ ' — ^ . (8) 



We show in Fig. ^ the total pressure on a layer of Al, 400 
nm in thickness. The pressure is calculated both with 
the stress-tensor method (solid grey line) and with the 
energy-gradient method (dotted grey line). The agree- 
ment between the two is excellent, demonstrating the va- 
lidity and usefulness of the energy-gradient picture. After 
all, the physical basis for forces in nature is the lowering 
of the energy, and electromagnetism is no exception to 
this. Thus the sharpness of the spatial variations of the 
stress tensor (or the energy) is a measure of how strong 
and localised forces are. Smooth, continuous variations 
of Taf3 imply forces that vanish when the volume element 
they act on shrinks to zero size. Surface forces, in con- 
trast, are localised at the interface between media; they 
are associated with abrupt, discontinuous changes in the 
stress tensor and can be explained physically as spatial 
discontinuities in the energy across the interface. 



B. Homogeneous system - evanescent waves 

The next level of complexity involves homogeneous sys- 
tems which are exposed to evanescent wavefields. In such 
fields the light wavevector k has got an imaginary com- 
ponent so that the dispersion relation is always satisfied: 



kk, 



(9) 



where e = e' + ie" is the dielectric function of the medium 
of propagation. 

Our system is a dielectric layer (denoted as medium 
II) of thickness d, separating two dielectric half-spaces I, 
III (see Fig. ||), with the planar interfaces lying in the xy 
plane. For simplicity we choose all media to have zero 
absorption. A beam of light of wavelength A is incident 
from medium I at an angle 9i and a wavevector ki = 
ap+z-fi . Upon penetration into region II, the wavevector 
may pick up an imaginary component along z: 



4 = /?' 



-7" ; k2 = (/?,0,72), 



when iS'^ < e'j 



(10) 

(/?, 0,^72), when/32>4. (11) 

By choosing e'2 < e[ and adjusting 61 (since (3 = ki cos 6*1) 
appropriately, we can create an evanescent wave in re- 
gion II with the desirable imaginary momentum 72 . Here 
we concentrate on the case of medium II being vacuum 
{e'2 = 1). We discuss total forces on media I and III, not 
interface forces. 

When only travelling waves exist in all three media, the 
two half-spaces repel each other. In this case, the photons 
can be thought of as tennis balls: as they scatter off the 
walls of I and II, they exert a positive pressure on each 
wall. 
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However, when k2 has a sufRciently large imaginary 
component along z, the force changes sign and becomes 
attractive. In this case the tennis-ball analogy breaks 
down. Our results indicate that at the electrostatic limit 
A ^ d, the total force on III per unit area varies as 



7Kl73p-7f) 
I71+73P 



(12) 



Thus the force becomes zero when I73I — 72'. Curiously, 
this is the same condition for the angles of refraction 
^2, ^3 to be equal, since 



tane*,;] = y— 



i = {l,2,3}. 



(13) 



When 72 < I73I the force is repelling, and Qi > 03, just 
as in the case of a travelling wave scattering off a denser 
medium. On the other hand, when 72 > I73I, we have 
^2 < ^3, which is reminiscent of a travelling wave scat- 
tering off a dense/rare interface; but at the same time 
the total force is attractive, a fact for which there is no 
analogue in waves of real k. 

The magnitude of the attractive force is sensitive to 
the ratio of the wavelength A to the layer thickness d 
over which the evanescent fields decay before entering 
medium III. We find that for \/d ^ 5 the force has al- 
ready dropped by more than three orders of magnitude 
compared to its value at the electrostatic limit A 3> d. In 
Fig. ^ we plot the normal pressure Pz versus cos 9i for 
three cases, (i) e'^ < £3, (ii) e'^ = 63 and (iii) e\ > £3. For 
(i) and (ii) the attractive part of Pz displays a smooth 
maximum, whereas for (iii) it peaks sharply. In all cases 
the highest value of Pz at the electrostatic limit is of the 
same order of magnitude as the incident radiation pres- 
sure. For (i) and (ii) the field in region III is that of a 
travelling wave for all angles and therefore I73 1 > 0. How- 
ever for (iii) the field in III is evanescent when 9i > 60, Oq 
being the angle for total internal reflection at an interface 
between media I and IIl[|. At 6*0 we have I73I — 0. From 
Eq. (^2|) it becomes clear why a vanishing 73 results into 
a sharp peak in Pz. For moderate absorption in III, I73I 
does not quite reach zero at 0q and the strength of the 
attraction is reduced although the peak is still clearly 
seen: for example, the maximum drops by 40 % when 
£3 ~ 0.1. Since the EM force is proportional to the light 
intensity Iq, we can adjust the latter so that the force 
induced by the evanescent (or the travelling) wavefield 
supersedes the dispersion forces between the half-spaces. 
We have indicated the strength of the Van der Waals at- 
traction in the figure, calculated in the macroscopic Lif- 
shitz theory with retardation and non-additivity effects 



'^It is interesting to note that a sharp feature in Pz is present 
always at 60, provided e'l > £3. Its position in the 6 axis is 
independent of the dielectric properties of layer II. 



considered. Owing to the behaviour of the Van der 
Waals force ||l^ , only large enough separations will allow 
the photo-induced force to dominate. 



C. Dielectric crystals 

We now turn our attention to a crystal of dielectric 
(GaP) spheres arranged in a SC lattice. GaP has an es- 
sentially real, constant dielectric function (modelled here 
as € = 8.9) over the frequency range of relevance to this 
work (0.07-0.9 eV). The actual absorption of GaP is too 
small to have any sizeable effect on our force calcula- 
tions, as we later show. We consider spheres with a ra- 
dius r = 200 — 400 nm, and choose a lattice spacing 
I = 900 nm, so the spherical surfaces are separated by 
D = 100 — 500 nm. We often consider the reference 
medium in which the spheres are embedded to be air 
i^med = l)j although, for reasons of experimental feasi- 
bility, it is really a liquid medium we have in mind. Such 
a medium does not alter qualitatively the effects we de- 
scribe provided its absorption is low enough. For exam- 
ple, we show that water is suitable. Most of our results 
are for a discretisation mesh 10 x 10 x 10. The crystal 
is infinite in the x and y directions, and has thickness N 
unit cells in the z direction. 

For the EM properties of this crystal we expect the 
following behaviour: 

(a) Propagation of EM waves in an "effective medium" 
at large wavelengths A ^ I where light senses the crystal 
as an homogeneous dielectric medium. 

(b) Multiple scattering effects: i) appearance of for- 
bidden ranges of frequencies as A '--^ ^ where the mul- 
tiply scattered waves interfere destructively and anni- 
hilate each other inside the crystal; ii) scattering by 
specific crystal planes into non-specular directions as 
A~^ ~ |k + G|/(27r), G being a reciprocal lattice vec- 
tor. 

(c) Resonant features at frequencies corresponding to the 
EM eigenmodes of each sphere. 

We observe these effects in Fig. ||, where the band 
structure (dispersion relation) is shown alongside trans- 
mittance and force for one layer of thickness, normal inci- 
dence. At low energies lu is proportional to the wavevec- 
tor kz, just as in free- wave propagation, the constant of 
proportionality being the wave velocity of a homogeneous 
effective medium replacing the crystal. For a system of 
unconnected spheres like ours, the dielectric function ee// 
of the effective medium is well approximated within the 
Garnett theory |13|, 



»(l + 2^) 
1-A 



2e„ 



(14) 



/ being the volume fraction of the spheres in the crystal. 
For our system Eq. ( p^ yields £e// = 1.76. We find that 
at low energies (or wavelengths A > 10-Z) the force on the 
crystal is well approximated by the analytical expression 
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for the force on a homogeneous film of the same thickness 
and dielectric function Ce//- 

At higher energies band gaps open up and the prop- 
agation of all EM modes is forbidden inside an infinite 
crystal: all incident light will be reflected back when the 
frequency lies within a band gap. It is impressive that 
even for one layer of thickness the major band gap cen- 
tred at around 0.7 eV causes a low transmission. Conse- 
quently the pressure on the layer approaches the upper 
natural limit 2Pi„c for frequencies inside this gap. The 
lowest band gap requires a thicker sample before it can 
clearly show its signature in transmission and force spec- 
tra, but otherwise it affects the system's behaviour in the 
same way. 

Another consequence of multiple scattering is that un- 
der certain conditions, the scattered fields can interfere 
constructively and produce outgoing waves in additional 
directions to those of the specularly reflected and directly 
transmitted beams. We call this scattering off-diagonal, 
since when it occurs the wavevector of the scattered field 
picks up a 2D reciprocal lattice vector component. 



k' = k 



(15) 



Most often the energetics of the problem are such that 
k'^ is imaginary and all the off-diagonal fields are evanes- 
cent. However at certain ranges of frequencies a real k'^ is 
allowed and the off-diagonally scattered beams are travel- 
ling waves. For a SC lattice the minimum energy for this 
to occur is when k|| lies at the edge of the first Brillouin 
zone along the TX direction. Considering an incident 
beam with k = (k^, 0, fc^), the condition becomes 



fcsin6' 



TT Co 



I 



I sin 6 



(16) 



Eq. ( |l6| ) gives the lowest frequency for off-diagonal scat- 
tering for a crystal of lattice constant I and a wave inci- 
dent at angle 9. 

As one imagines, the effects on forces are rich as the 
possible combinations of momentum-carrying outgoing 
waves are endless. To give one example, consider the 
incidence at 9 = 70deg of a laser beam on our crys- 
tal of one layer thickness. For this system we plot the 
intensity of reflection r and transmission t of each out- 
going wave in Fig. |6[ Off-diagonal scattering sets in for 
511 — {~2n/l, 0, 0) at frequency 0.71 eV in good agree- 
ment with Eq. (|l6|). At slightly higher frequencies the 
extra energy is given to k'^ while k'^ is fixed to the value 
tt/1, provided the angle of incidence is not changed. Thus 
the outgoing waves include two reflected and two trans- 
mitted beams, all forming an angle 9 with the surface 
normal (Fig. In addition, notice that at 0.76 eV the 
intensities of the leftwards and the rightwards travelling 
waves are equal; + = + . This means that 
the total momentum of the scattered EM fields along x 
is zero. In other words, the crystal absorbs all incident 
momentum along x, and the a;-force on a surface S equals 



= iP^ncS) sm9. 



(17) 



Our force calculations confirm that approaches this 
value at w ~ 0.76 eV. 

Off-diagonal scattering survives a moderate amount of 
absorption in the crystal. We find that for e" = 0.1 
the reflectivities of the off-diagonal beams decrease by 
less than 40% (depending on the frequency), which is a 
little larger than the decrease the diagonal beams suffer 
(< 20%). 

Thirdly, we note that at specific frequencies there ap- 
pear highly flat modes in the band structure combined 
with sharp features in the transmittance spectrum (for 
example at 0.624,0.768,0.776,0.786,0.85,0.88 eV on the 
kz = axis. Fig. Due to their low dispersion they 
possess a low group velocity duj/dk and for this rea- 
son they have been conveniently termed 'heavy photons' 
by Ohtaka and Tanabe |14|. As has been argued else- 
where ^4 15 1 , these are resonance peaks corresponding to 
the EM eigenmodes of isolated spheres {Mie resonances) . 
Their presence strongly affects the position of the band 
gaps on the frequency axis |l^. The eigenmodes are 
spherical waves labelled by three quantum numbers, just 
as in the case of electronic orbitals in atoms, and are gen- 
erally electric or magnetic oscillating multipoles. Due to 
the spherical symmetry they possess a (2n + 1) degener- 
acy, n being the azimuthal quantum number. However 
the presence of a plane wave incident along the direction 
of k breaks this symmetry and modifies the boundary 
conditions so that only one eigenmode for each n can be 
excited; this, at any rate, is the standard Mie theory re- 
sult and strictly applies to the case of an isolated sphere. 
The notation for the excitable electric (magnetic) modes 
is (bli), i being the sequential index of the resonance 
in the frequency axis. For dielectric spheres the electric 
modes a\ are broader whereas the magnetic modes are 
sharper, the sharpness of both types of modes increasing 
with frequency. For metallic spheres, only the electric 
modes survive, being better known to some as 'surface 
plasmons'. The first few Mie resonances for one GaP 
sphere subjected to a plane wave are shown in Table |. 

When many spheres approach each other to form the 
lattice, interactions occur which shift the resonances with 
respect to the frequencies for an isolated sphere (Fig. 

Unless the coupling between neighbouring spheres 
is strong, each resonance in the crystal will retain its 
characteristic EM field distributions on the surface and 
the interior of each sphere. Thus the lowest-in-frequency 
heavy-photon resonance in our GaP crystal is the h\ mag- 
netic dipolc mode, followed by the electric dipole a\ mode 
(Fig. 1^). The EM fields for these two resonances clearly 
match the expected patterns for Mie scattering from a 
single sphere. For the next two sharp peaks in trans- 
mittance the EM fields are similar but not identical to 
those for the magnetic quadrupole mode h\] and likewise 
for the following two peaks, whereby the fields display 
partial similarity with those for the electric quadrupole 
mode a\. These are cases of mixing from more than one 
Mie modes owing to the organisation of the spheres into 
a lattice. To quantify the contribution of each Mie mode 
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into the heavy-photon resonance one has to decompose 
the scattered field into its constitutive spherical waves, as 
Ohtaka and Tanabe have done [|^. However the mode 
that follows in frequency is a distinct electric quadrupole, 
with field distributions identical to the 02 Mie mode. 

Depending on lo, the incident field produces polarisa- 
tion charges and/or induced currents in the spheres and 
orchestrates their response according to which mode is 
excited each time. It is interesting that spheres which 
are initially non-magnetic (/x = 1) can acquire magnetic 
properties at the occurrence of any &^ resonance. An ef- 
fective homogeneous medium replacing the crystal struc- 
ture would have /i 7^ 1 at this instance. 

What happens at a resonance? Each sphere acts as 
a cavity that traps the EM waves. The energy density 
internal to the spheres is 2-3 orders of magnitude larger 
than for the incident light beam As Van Albada 

et al have shown the transport of EM energy can 
be slowed down considerably at the conditions for Mie 
scattering. They derive a rigorous microscopic formula 
for the velocity of energy transport ve, and subsequently 
propose a heuristic but physically plausible model for ve, 
which they show to agree with the exact model and ex- 
perimental results within reasonable accuracy, and which 
reads 



Ve 



Co 



(18) 



where cq the speed of light in vacuum, / the volume frac- 
tion and W the energy density internal to the spheres rel- 
ative to the surrounding medium (vacuum in our case). 
Our computations suggest that ~ 20 for the first few 
heavy-photon modes, resulting in a velocity ve that is 
much lower than one expects from a straightforward es- 
timate for the average velocity away from resonances, 



Co 



/(v^-l) + i 



(19) 



Thus for / = 0.22 the equations above give ve/cq ^ 0.19, 
whereas Vav/co ~ 0.69. For comparison, the group veloc- 
ity Vg — duj/dk varies with frequency and ranges from 
0.02 Co to 0.25 Co for these resonances. The accumula- 
tion of energy in the system at the resonance condition 
means that the overall EM fields are enhanced^, and this 
results in large forces when the stress-tensor integration 
is carried out. 

The magnitude of the photo-induced force is propor- 
tional to /q, the intensity of the incident beam. Pro- 
vided the spheres are in a liquid medium which facil- 
itates heat evacuation (so that melting is avoided), the 



photo-induced force can dominate gravitational, thermal- 
fluctuation and Van Der Waals forces. Assuming a mod- 
erate laser of Jo '--^ 3.5 x 10® W/m^, the incident radiation 
pressure is ~ 1 pN per unit cell (900 x 900 nm^). For a 
crystal of one layer thickness, and when lu lies inside the 
second band gap (centred at 0.7 eV in Fig. ^) the force 
on a unit cell is ~ 2 pN (the 'natural limit' for full re- 
flection), whereas at resonance conditions it exceeds this 
limit and can reach much higher values (e.g. 5 pN for the 
b\ resonance). By contrast, the gravitational force on a 
GaP sphere is mg ~ 0.8 x 10~^ pN. Thermal-fluctuation 
forces are comparable to gravity. Finally, the Van der 
Waals attraction between two GaP spheres^ is less than 
0.2 pN |l|. 

The sharp force features shift slightly with the angle 
of incidence. More importantly, they are sensitive to the 
amount of absorption present in the spheres or the sur- 
rounding medium: a moderate e" = 0.1 is enough to 
suppress them. For GaP spheres in water, however, we 
have found that the absorption is too small to blur the 
sharp features in forces, within the frequency range of 
interest (Fig. |l|). 

The effect of increasing thickness is of particular impor- 
tance with reference to resonances. When the number of 
layers is doubled, the number of sharp features doubles 
also [|l^. This effect is directly equivalent to degener- 
acy splitting in atoms. Two oscillating multipoles (GaP 
spheres in a heavy photon mode) brought close together 
interact, their fields hybridise and the two new energy 
levels of the combined system lie slightly below (at uj^) 
and slightly above (at uj^) the frequency of the original 
level (at uq). This atomic physics analogy is far reaching: 
the mode at lu^ is bonding (attractive forces) whereas at 
a;> it is anti-bonding (repulsive forces). We see a series 
of alternating signs in the force on a sphere as different 
resonances are encountered (Fig. |ll|) for a crystal two 
layers thick. For the magnetic dipole {b\) mode in par- 
ticular, we have shown Pq ] that at a;< all spheres have 
their dipoles parallel, whereas at w> the dipoles of the top 
layer are anti-parallel to those of the bottom layer. The 
same is true for the electric dipole (aj) mode: all dipoles 
are parallel at a; = 0.7915 eV, anti-parallel a.t uj — 0.8025 
for the system of Fig. |ll|. Similarly, the multipoles of 
the higher-energy modes can be oriented constructively, 
in which case an attraction occurs, or destructively, lead- 
ing to repulsion between neighbouring spheres. This ar- 
gument holds also for the mixed hybrid states that are 
made out of more than one Mie modes, for example at 
0.86 eV and 0.875 eV in Fig. |ll|. 

As expected, the splitting of the energy levels and the 
force strength increase as the spheres of two vertically 
adjacent layers approach each other (Fig. |l^). The na- 



■^This mechanism has been known to explain the enhance- 
ment of the field on the surface of GaP particles used in SERS 
(Surface Enhanced Raman Scattering) experiments |2c|]. 



■^This value is obtained in the D <^ r limit and ignoring 
retardation effects and is therefore an overestimate. 
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ture of the sharp forces between adjacent layers therefore 
Ues in the sphere-sphere interactions induced by the hght 
beam, when the latter populates the various Mie reso- 
nances. 

The analogy between resonant, light-induced, sphere- 
sphere interactions and the bonding of atoms that we just 
described is not the only one to be made. An analogy 
exists also with the case of cutting a metal in two along a 
plane, whereby the bulk plasmons reorganise into surface 
plasmons of attractive and repelling nature at frequencies 
below and above the original frequency The similar- 
ities between these three very different systems prompts 
us to acknowledge a general tendency in nature for the 
formation of energetically favourable EM bonds between 
two systems in the same state that come into proximity. 

The rest of the force spectrum displays mostly smooth, 
attractive interactions due to the band gaps between the 
heavy-photon bands. This off-resonance attraction is 
produced as the band-gap condition for high reflectivity 
causes the top layer to absorb almost all the shock from 
the reflection of the incident light and be therefore pushed 
towards the bottom layer which remains relatively inert. 

Next we explore the possibility that the induced bond- 
ing interactions may lead to the formation of stable crys- 
talline structures with full 3D periodicity. We start with 
two crystal layers on top of each other in a SC arrange- 
ment and gradually displace the top layer along the sur- 
face diagonal (1,1,0), such that the relative distances 
between spheres of the same layer do not change, as 
shown in Fig. |l3|. We monitor the normal (Fz) and shear 
{Fx,Fy) forces for the bonding mode of the magnetic- 
dipole b\ resonance. 

We find that relative to a given sphere of the bottom 
layer (^4), the shear force on the sphere above {Ai) al- 
ternates in sign as the top layer is displaced along the 
diagonal. For small displacement S < 0.3, the sphere 
Ai is attracted back to its initial position Sinit — 
i.e. directly above A of the bottom layer. As S in- 
creases (0.3 < S < 0.7) the top sphere is attracted to- 
wards the void region in the centre of the square ABCD 
{Smid = 0.5), with a shear- force strength 5 times as large 
as the attraction to Sinit- Finally for 0.7 < S < 1.0 it 
is attracted towards Sfm — 1.0, i.e. directly above the 
sphere C of the bottom layer, with a force strength sim- 
ilar to that for the attraction towards Sinit- The normal 
force between the two layers remains attractive (since 
we are monitoring the bonding mode) but its magnitude 
displays a global maximum at Smid and two local max- 
ima at Sinit ) S f in . Consequently, the spheres of the top 
layer experience two potential wells at Sinit , S f in , whereby 
a stacking sequence AAA... is obtained. However these 
wells are relatively shallow and when the potential bar- 
rier is crossed, the spheres find themselves in a much 
deeper well at Smid where a sequence ABAB... is ener- 
getically favourable (Fig. |l^). As for the anti-bonding 
mode, there are two high potential barriers at Sinit and 
S fin and a barrier of lower height at Smid - 

The occurrence of resonant bonding effects between 



spheres as well as the energetic preference for the arrange- 
ment of crystal layers into the sequence ABAB..., may 
add to the current methods for engineering 3D photonic 
crystals out of assemblies of colloidal spheres in liquid so- 
lutions. Experiments on gold nanospheres (radius r ~ 10 
nm) by Kimura have shown fast formation of aggregates 
upon illumination at the surface-plasmon frequency [i.e. 
the electric-dipole mode aj), in agreement with calcula- 
tions of the photo- induced Van der Waals dipole-dipole 
interaction p3]. Burns, Fournier and Golovchenko on the 
other hand have concentrated on dielectric microspheres 
(r ~ 700 nm) p3t . They have experimentally demon- 
strated the presence of long-range attractive forces, in- 
duced by laser light, between two polystyrene spheres in 
water. They argue that these are magnetic-dipole in- 
teractions, owing their existence to induced currents on 
the spheres. They calculate the interaction energy be- 
tween two magnetic dipoles based on a simplified model 
that works in the long-distance limit R ^ X, where R 
the distance between the dipoles and A the wavelength 
of light. From the perspective of Mie resonances this ap- 
proach is justified because although the spheres Burns 
et al use possess many different resonances at or near 
the frequency of their laser light, the large inter-sphere 
distance ensures that the dipole approximation works to 
leading order. Also the electric-dipole solutions are of 
secondary importance since they are damped compared 
to the magnetic dipoles. 

In contrast to the work of Kimura (metallic spheres) 
and that of Burns et al [R ^ A), we are dealing with 
dielectric spheres at close proximity {R ^ A). In gen- 
eral, the sphere-sphere interactions in our system display 
a mixing of EM-multipolc behaviour, unless a specific 
condition for an isolated Mie resonance is met. For ex- 
ample two spheres of r = 349 nm interact at ^ 0.65 eV as 
oscillating magnetic dipoles, at ^ 0.80 eV as oscillating 
electric dipoles and so on. Although the model of Burns, 
Fournier and Golovchenko does not apply to our system 
(due to the close proximity between the spheres) , we be- 
lieve that their physical explanation for the spatially os- 
cillatory shear forces is basically correct: the oscillation 
is due to the changes in the phase shifts associated with 
retardation between spheres of different layers. 

IV. CONCLUSION 

In conclusion, we have presented a general methodol- 
ogy for accurate and efficient computation of EM fields 
and forces in matter. We have applied this first to ho- 
mogeneous systems (dielectric and metallic) and stud- 
ied the effects of travelling and evanescent waves on the 
forces exerted on a slab of material. Under specific condi- 
tions, evanescent fields can produce attractive total forces 
which are even stronger than the Van der Waals attrac- 
tion between two half-spaces. Next, we have concen- 
trated on a dielectric crystal of GaP spheres arranged in a 



8 



SC lattice. We observe the connection between lattice ef- 
fects (band gaps, Bragg scattering) and EM forces. More 
interestingly, however, we find large resonant forces at 
frequencies corresponding to the EM eigenmodes of iso- 
lated spheres. An analogy to atomic physics proves fruit- 
ful as the EM fields on the spheres arrange into bonding 
and anti-bonding states, causing highly attractive and 
repelling interactions which may become dominant over 
gravitational, thermal and Van der Waals forces. We be- 
lieve these effects to be experimentally observable and 
relevant to the formation of stable crystals of colloidal 
nanospheres. 
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FIG. 1. The intensity of reflected ligfit, at normal inci- 
dence, from two semi-infinite fiomogeneous systems, GaP and 
Al, is shown in black lines. The sharp decline in reflectance 
for Al is due to the plasma frequency ujp = 15 eV. A fi- 
nite but substantially thick (400 nm) layer of Al also reflects 
most light below ujp; above ujp it displays a ripple structure 
due to multiple refiections of the waves from the sides of the 
layer (Fabry-Perot oscillations). Reflectance (transmittance) 
for the finite layer is shown in solid (broken) grey lines. 



FIG. 2. Pressure on a homogeneous medium of infinite 
(GaP and Al half-spaces, black curves) and finite (Al 400 
nm, grey curves) thickness. Light is normally incident along 
the +z direction with intensity Jo ~ 3 x 10* W/m'^. Positive 
values for the pressure imply that the light is pushing the sys- 
tem, (i) For the half-spaces, the pressure is calculated with 
the stress-tensor method analytically (solid black lines) and 
numerically (dotted black lines) in perfect agreement. It lies 
within Pine (full absorption) and 2Pi„c (full reflection) irre- 
spective of the different dielectric properties of the GaP and 
Al media, (ii) In contrast, a layer of finite thickness experi- 
ences a pressure with a zero lower bound, corresponding to 
full transmission through the layer. For this system the pres- 
sure is calculated numerically with the stress-tensor method 
(solid grey line) and with the energy-gradient method (dotted 
grey line). 



FIG. 5. (a,b) Photonic band structure along the FX di- 
rection for a simple cubic lattice of GaP spheres (r — 365 nm, 
e = 8.9) in air. The vertical axis is in eV units for the fre- 
quency, (c) Transmittance for a crystal one-unit-cell thick and 
for normal incidence along z. The sharp peaks correspond in 
frequency to the low-dispersion modes of the band structure 
in (b), and occur when the incident light excites one of the 
EM eigenmodes of isolated spheres (Mie resonances). Some 
peaks can be attributed to a single resonance (6i: magnetic 
dipole mode; a\: electric dipole mode; electric quadrupole 
mode) while others involve mixing from more than one Mie 
modes (hybrids) owing to the organisation of the spheres in 
a lattice, (d) Normal force on one unit cell for the system of 
(c) and beam intensity /o = 3 x 10* W/m^. The force scales 
with Jo. 

FIG. 6. Intensity of the specularly refiected (solid black 
line), directly transmitted (solid grey line), off-diagonally 
refiected (dashed black line) and off-diagonally transmit- 
ted (dashed grey line) beams. The angle of incidence is 
9 — 70 deg. For low frequencies, corresponding to wavevectors 
away from the edge of Brillouin zone, only the diagonal beams 
exist. However at the condition of Eq. $lm off-diagonal scat- 
tering starts to occur (at 0.71 eV). For the frequency range 
shown here only two off-diagonal beams exist. The spheres 
(e = 8.9) have radius 365 nm and are suspended in air. 



FIG. 3. A wave incident at angle Oi from medium I is scat- 
tered by the interfaces with media II and III, which are sepa- 
rated by a thickness d. When 62 < ei the fields in II become 
evanescent for large enough 9i. 

FIG. 4. Total pressure on medium III (normal force per 
unit area), plotted against the angle of incidence, for (i) 
e'l = 8,e^ = 1,4 = 9, (ii) e[ = 9, = 1, e?, = 9 and (iii) 
e'l — 10,62 ~ 1,£3 ~ 9. For the latter system nonzero ab- 
sorption (£3 = 0.1) is considered also (grey line). Negative 
pressure implies attractive force on III. The sharp attraction 
peak for (iii) occurs when I73I = 0. The zero pressure corre- 
sponds to 72 = I73I. All curves are for d — 100 nm, X/d = 100 
and an incident beam intensity Jo — 3.5 x 10^ W/m^. For 
this system the Van der Waals attraction is smaller than the 
photo-induced force for most angles; however for light inten- 
sities weaker by a factor of 20 or more, or for separations 
smaller than 50 nm the Van der Waals force dominates all 
photo-induced effects over all angles. 



FIG. 7. Geometry of the scattered beams for off-diagonal 
scattering. The off-diagonal beams have negative momentum 
along X. At 0.76 eV, when the intensities of the waves travel- 
ling to the right and to the left balance each other (see Fig. |^, 
the total scattered a;-momentum is zero; the crystal absorbs 
all incident momentum along x. 

FIG. 8. For a given Mie mode (at frequency lu) of an iso- 
lated sphere (radius r), the order parametre u ■ r is indepen- 
dent of the sphere size. When the spheres are put together in 
a lattice, interactions shift the resonances. When the radius 
is small compared to the separation {i.e. the unit-cell length 
I), the one-sphere result is recovered. Results here are for the 
b\ mode. 

FIG. 9. (a): The spherical sections Si, 82,83 on which 
the field distributions are plotted, (b-d): Dipole mode. For 
the bl resonance (magnetic mode) the H field is plotted on 
the spherical sections 81,82, while the E field is shown on 
83. For the a\ electric mode the role of the H and E fields is 
interchanged, (e-g): Quadrupole mode. For the hi resonance 
(magnetic mode) the H field is plotted on the spherical section 
Si, while the E field is shown on 82,83. For the 02 electric 
mode the role of the H and E fields is interchanged. 
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FIG. 10. Dependence of the force on absorption. Away 
from resonances the forces do not suffer a sizeable change 
with moderate absorption. The resonances are sensitive to 
losses, but the actual amount of absorption present in GaP 
spheres and in an experimentally suitable liquid medium (in 
this case water) is small enough so the resonance effects can 
be clearly seen. Here we plot the normal pressure on a crystal 
one layer thick (GaP spheres, radius 337 nm) in air (dashed 
grey line) and in water (dashed black line) when a light beam, 
of intensity 7o = 3 x 10* W/m^, is incident along z. For 
comparison, the force on a hypothetical loss-free GaP sphere 
in air is also shown (solid black line); it is clear that GaP is 
a material with very low absorption in this frequency range 
and is therefore highly suitable for displaying the effects of 
resonant scattering. 

FIG. 11. Normal pressure on a two-layer crystal of GaP 
spheres (e = 8.9, radius r = 349 nm) in air when light is 
incident along z with intensity 7o = 3 x 10* W/m^. The 
pressure on the top (bottom) layer is shown in solid black 
(grey) line; on both layers {i.e. total pressure) it is shown in 
a dotted line. At resonant frequencies the total pressure dis- 
plays peaks of mostly positive (repelling) pressure. However 
the resonant forces on either layer are generally opposite in 
direction, peak more sharply than the total force and their 
sign alternates with frequency, exhibiting a series of bond- 
ing/anti-bonding interactions (i.e. layers attract/repel each 
other). At the bonding modes the spheres of the two layers 
behave as oscillating EM multipoles oriented parallel to each 
other, at the anti-bonding modes they are anti-parallel. Vari- 
ous modes are seen: magnetic dipole (&}), electric dipole {a\), 
electric quadrupole (al) and hybrids. Compare with the force 
on one layer of Fig. ra. 



FIG. 14. The potential landscape experienced by sphere Ai 
while it and all other top-layer spheres are being displaced 
along the surface diagonal AC relative to the bottom-layer 
spheres. The potential wells at Si„it = and = 1.0 

lead to a stacking order AAA..., whereas the deepest well at 
Smid ~ 0.5 favours a sequence ABAB.... The potential energy 
is in units of 10~^* J. Results here are for the bonding b\ mode 
of spheres of radius r = 240 nm, e = 8.9 and a lattice constant 
/ = 900 nm. 

TABLE I. Mie modes of an isolated GaP sphere (e = 8.9) 
and their corresponding order parametre (mode frequency 
times sphere radius) in units of eV x nm. (aj,) modes rep- 
resent magnetic (electric) multipoles oscillating at the light 
frequency. 
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FIG. 12. Resonant forces on a two-layer sample when light 
is incident normally along z with intensity Jo = 3 x 10* W/m^ 
and populates the magnetic dipole mode b\ . As the two crys- 
tal layers approach each other the bonding/anti-bonding en- 
ergy levels of the Mie modes on the spheres split further, and 
the potential well (barrier) for attraction (repulsion) between 
the layers becomes deeper (higher). The pressure on the top 
(bottom) layer is shown in black (grey). Attraction occurs 
when the pressure on the top layer is positive and that on the 
bottom layer is negative. The numbers denote the distance 
between two vertically adjacent spheres in multiples of the 
unit-cell length, I = 900 nm. The GaP spheres have radius 
r — 365 nm and are in water. The absorption of both the 
spheres and the medium is taken into account. 



FIG. 13. View from above of a sample two layers thick. All 
the spheres of the top layer (white) are displaced along the 
diagonal AC while the bottom-layer spheres (grey) remain 
fixed. The position of a top-layer sphere can be traced by the 
variable 5 = x/l, I being the lattice constant for the original 
simple cubic lattice. 
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